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Motivated by recently developed techniques making it possible to compute Casimir energies for 
any object whose scattering S-matrix (or, equivalently, T-matrix) is available, we develop a variable 
phase method to compute the S'-matrix for localized but asymmetric sources. Starting from the case 
of scalar potential scattering, we develop a combined inward/outward integration algorithm that is 
numerically efficient and extends robustly to imaginary wave number. We then extend these results 
to electromagnetic scattering from a position-dependent dielectric. This case requires additional 
modifications to disentangle the transverse and longitudinal modes. 

PACS numbers: 03.65.Nk, 11.80.Et, ll.80.Gw 



Scattering theory Jlj, [2( is an invaluable tool for investigating a wide range of physical systems. Far away from a 
system that is localized in space, one can express solutions to the wave equation as free incoming and outgoing partial 
waves. In this partial wave basis, the scattering S'-matrix then gives the amplitude and phase of outgoing waves 
reflected from the system in terms of a given amplitude and phase of incoming waves. 

One of the many applications of scattering theory arises in calculating Casimir forces. While the connection 
between Casimir forces and scattering amplitudes has long been understood in planar systems 0, 0| , only recently 
have techniques been developed in which the Casimir force is expressed in terms of the S'-matrix (or, equivalently, T- 
matrix) for general geometries 043] • In this approach, the S-matrix encodes the effects of quantum fluctuations on a 
single object, while universal translation matrices, obtained from the free Green's function, encode the objects' relative 
positions and orientations. This decomposition provides a concrete implementation of the "TGTG" representation of 
the Casimir energy in terms of scattering transition operators and free Green's functions Q . The S-matrix is also a 
key ingredient in Casimir calculations of quantum corrections to soliton energies and charges These calculations 
take advantage of the relationship between the S-matrix and the change in the continuum density of states, 



where the eigenvalues of the matrix in parentheses are the scattering phase shifts. 

A standard approach to finding the exact electromagnetic S-matrix for dielectric objects involves integrating the 
vector solutions of the Helmholtz equation in dielectric media over the object's surface [Ijj, LLLl ■ A variety of subsequent 
techniques have obtained a wide range of analytic and numerical results 0. In myiases of practical interest, 
one can obtain approximate results valid in appropriate limits, such as large or small values of the wavelength or 
partial wave number. Because the Casimir calculation involves summing over all fluctuating modes, however, suitable 
approximations are often not available. 

In Casimir problems and other applications of scattering theory, one frequently considers objects with sufficient 
symmetry that the problem separates and the S-matrix is diagonal. For cases where the resulting ordinary differential 
equation cannot be solved analytically, the variable phase method [1, [l]| provides an efficient numerical algorithm. 
In particular, it allows one to solve for the S-matrix as an initial value ODE, rather than a boundary value problem. 
Here we extend the variable phase method to compute the S-matrix in situations without any symmetry assumptions. 
We begin with the case of a scalar potential, as arises in quantum mechanical potential scattering, which we then 
generalize to the case of electromagnctism with a position-dependent dielectric. Our approach can provide a middle 
ground between anal ytic results and fully general numerical calculations For dielectrics, our work extends the 
results found in Ref. [l5j], which treats the case of a spherically symmetric but r-dependent dielectric. As we will see, 
the asymmetric case introduces additional complications, because one can no longer rely on the channel decomposition 
to separate transverse and longitudinal modes. We also introduce a combined inward/outward integration algorithm, 
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which makes use of the Wronskian of the regular and outgoing solutions, to ensure the stability of the numerical 
calculation for imaginary wave number k = in. 



II. HELMHOLTZ SCATTERING 



We begin by considering scattering of waves obeying the scalar Helmholtz equation, as would arise in a typical 
quantum mechanics problem. This calculation generalizes straightforwardly to the vector Helmholtz equation, as we 
show in this section. Additional formalism is needed for the case of Maxwell scattering, however, so we postpone that 
case to the next section. 



A. Variable Phase Approach: Outgoing Wave 



We start from the Helmholtz equation in three dimensions 

- V 2 4'k{r) + V(r)Mr) = k 2 ^ k {r) , (2) 

where the potential V(r) is localized in a region around the origin. This equation describes, for example, ordinary 
quantum-mechanical scattering of the scalar wavcfunction ipk(r) from a localized potential. Since each k value is 
treated separately, V(r) can also be fc-dependent, though we do not indicate this possibility explicitly. We expand 
both the solution i/jk{f) and the potential V(r) using a Fourier series in the angular variables, 

oo i 



*«=EEjw r ) F ™(^) and v ( r ) = J2 E w(r)y/(M), (3) 

i=0m=-( e'=Om'=-i' 

to obtain 

]T Y e m (9, 4>)(~j£s + - k 2 ) ^m,k{r) + lfcm,*(r)Y?*(0, 0) £ V t . m >{r)Y t ?'(6, 0) = . (4) 

im ^ ' tm I'm' 

Next, we multiply both sides by YgJ} (0,0)* = (— l) m Y^, m (0,0) and integrate over solid angle. The last term 
becomes a convolution, which mixes angular momentum channels. We obtain 

+ l"{t" + 1) _ fc2 \ ^ /fm/fifc(r) + ( £ V t , m ,(r)Z%#, m " J W (r) = , (5) 



Km \t m 



where the integral identity 



2 \, vra/fl .xW//, W /fl ,s /(2^+l)(2£' + l)(2f + 1) (I £' l"\ ( £ £' t" . 

#yi m (0,0)y;i™ (0,0)^ (0,0) = 4?r ( ) [ m m , mtt ) (6) 
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allows us to express Z$™„ m in terms of 3j-symbols as 



Z "'™' m _ ( 1)m V 4^ [OO ) [m m! ~m") ' (7) 

In the absence of a potential, the regular and outgoing solutions for V&n,fc( r ) are given in terms of spherical Bessel 

and spherical Hankel functions by krjg(kr) and krh/fc (kr), respectively. 

Since the scattering channels will mix for a nonspherical potential, we will want to consider all incoming waves 
together. To do so, we rewrite Eq. (|S|) as a matrix differential equation, 

(-^5 + ^ - fc2 ) Mr) + V(r)Mr) = 0, (8) 

where hat indicates a matrix indexed by the angular momentum indices £ and m (so that both i and m are combined 
into a single matrix index), L 2 is a diagonal matrix with £(£+1) on the diagonal, and V(r) is the matrix in parentheses 
in the second term of Eq. ([S]). 
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We begin by considering the solution to this equation with outgoing wave boundary conditions, which we parame- 
terize as 

F k (r) = G k {r)W(kr) , (9) 

where W(x) is a diagonal matrix with the free outgoing wave solutions xhf\x) on the diagonal. The Helmholtz 
equation for F k (r) then translates into an ordinary differential equation for the matrix G k (r), 

-GUr)-2Gi(r)(^\ogW(kr))+\[L 2 ,G k (r)}+V(r)G k (r)=0, (10) 



\dr 

where prime denotes derivative with respect to r and we have used the fact that the free solution obeys 

L 2 

- W"{kr) + — W(kr) = k 2 W{kr) , (11) 

and then multiplied from the right by W (kr). By the outgoing wave boundary condition, we have G k (oo) = 1 and 
G' k (oo) = 0, where 1 and are the identity and zero matrices respectively. These results provide the necessary initial 
values for integrating Eq. (|10[) inward from infinity to the origin. 

To define the S- matrix, we combine the solutions with k and —k (or, equivalently, the outgoing wave solution and 
its conjugate, the incoming wave solution) to form the physical wave function, 

Mr) = -G. k (r)W(-kr)M + G k (r)W(kr)S k (k) , (12) 

where M is a diagonal matrix with (— l) e on the diagonal. We then find the S-matrix by the regularity condition at 
the origin, which yields 

S k = lim#- 1 (fcr)Gr 1 (r)G_ fe (r)VK(-fcr)M. (13) 

r— >0 

In many applications it is convenient to work with the T-matrix, which is given by T k = ^(S k — 1). 

We can thus find the S'-matrix numerically, by integrating G k (r) in from r = oo to r = 0, and similarly for G— k (r). 

d 

The combination W' (x)W~ 1 (x) = — — log W(x) is easy to calculate numerically, since it is just a diagonal matrix with 

ox _ _ . — . 

rational functions of x on the diagonal, which can be obtained from a finite continued fraction expansion |16l . p. 241]. 

The inputs to the calculation are then the "multipole moments" of the potential at each r, V( m (r). We could imagine 

some simple non-spherical potentials for which these moments might be particularly easy to find; or, we could specify 

the potential explicitly through its representation in this spherical harmonic basis. 

Note that in the ordinary variable phase method, where the channels separate (so here all matrices would be 

diagonal), it is common to write G k (r) = e l/3k ^ r \ which further simplifies the calculation. This approach is problematic 

in the general case, however, because then fi k (r) doesn't commute with its derivatives. 



B. Variable Phase Approach: Regular Wave 



In principle, one could carry out the calculation of the previous subsection for k = in to obtain the S'-matrix on the 
imaginary fc-axis, as is typically required in Casimir calculations. In practice, however, this is not possible, because in 
place of the oscillating spherical Bessel function h^\kr), we now have the exponentially decaying modified function 
ke(Kr), which then grows exponentially as we integrate in from infinity. As a result, a direct application of the previous 
results is hopelessly unstable numerically, and we will need to introduce some additional formalism to obtain a useful 
calculation. 

To address this problem, we use an approach developed in Ref. (l7| , in which we parameterize the regular solution in 
a complementary way to what we did for the outgoing solution in Eq. ^ . Here it will be convenient to parameterize 
the transpose of the regular solution, as 



h(rY = W{kr)- x H k {r) 
(note the reversed order in this decomposition). We then have 



^logV^(fcr) 



- [L 2 ,H k (r)] + H k (r)V(r) = 0. 



(14) 



(15) 
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By the regularity of $fc(r)* at the origin, we have the boundary condition H k (0) = and H' k (0) = 1, where again 
prime denotes a derivative with respect to r. Starting from this boundary condition, we can then integrate Eq. (|15[) 
outward from the origin. 

This integration also contains instabilities for k imaginary, but what will be useful to us is that they show up in 
a complementary region: The integration of G k (r) blows up for r — > 0, while the integration of H k (r) blows up for 
r — > oo. We can make use of this complementarity by considering the Wronskian of our two solutions [l|, p. 465], 

mi = W[$ fc (r)*,&(r)] = Mrf - (J^M*) A(r) 



= [W(kr)- 1 H k {r)j [G k {r)W \kr) + G' k {r)W(kr) 

3 (W(kr)- 1 ^ H k (r) +W(kr)- 1 H k {r)^j (& k (r)W(kr) 

H k {r) (d k (r) (-?- \og W(kr)) + G' k (r) 



dr 
Wi^kr)- 1 







H' k {r)~[—\ogW(kr))H k {r))G k {r) 



W(kr) , 



which is independent of r. By the boundary conditions on G k {r) and H k (r), we also have 

limW[$fe(r-)*,F fc (r)] = lim \-W{kr)- l G k {r)W{kr) 

Thus, at any r, 



r->0 



W[$ k (ry,F k (r)} 



lim 

r-s-0 



-V^(fcr)- 1 G fc (r)W^(fcr) 



(16) 



(17) 



(18) 



But the right-hand side of this equation gives the quantity we need to calculate the S- matrix from Eq. (|13p . So our 
strategy will be to pick an intermediate radius r$ and integrate both G k (r) in from r = co to r = ro and H k (r) out 
from r = to r = r^. Then we can evaluate the Wronskian in Eq. (| at r = ro and use it to obtain the right-hand 
side of Eq. (|18|) , which is what we need to find the S'-matrix. This procedure will continue to be stable even when k 
is imaginary (with either sign of its imaginary part — and we will need both signs to compute the S- matrix) . 
We thus obtain 



§k=(w[*k(r) t ,F k (r)]- 1 (W(kr)- 1 W(-kr) (w[$_ fc (r)*,F_ fc (r)] 
This expression is now suitable for numerical evaluation. 



M. 



(19) 



C. Vector Helmholtz Equation 



We next generalize this calculation to the vector Helmholtz equation, 

- VVfc(r) + V(r)xp k (r) = k 2 ^ k {r) , 



(20) 



where our wavefunction is now a three-component vector ip k {r). Our eventual goal is to study electromagnetic 
scattering, which will require significant additional modifications of this approach to disentangle the transverse and 
longitudinal modes. In contrast, the generalization to the vector Helmholtz equation is relatively straightforward, 
requiring only that we establish corresponding definitions and identities appropriate to the vector case, which we take 
from Ref. [l|. 

We begin by defining the three vector spherical harmonics for each value of j = 0, 1, 2, 3 . . . and m = —j . . . j, 



+i 



(21) 



a— — 1 m' — - 



where I = j, j ± 1 for our three vector spherical harmonics, Cf™ la is a Clcbsch-Gordan coefficient, and the spherical 
basis vectors are 



ei = -= I sin r + cos + i<j> 

v2 ^ 



1 

71 



x + iy) 
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eo = cos 6 f — sin 6 = z 

e _i * / ~\ 1 

e_i = [sm6r + cos96 - i(f>j = — (x - iy) . (22) 

For j = 0, we have only the case £ = 1. This representation effectively couples the orbital angular momentum £ to 
the s = 1 spin angular momentum associated with the vector index. We can then decompose tp{r) as 

oo j+1 j 

^W = E E E -*iM«(M). (23) 

j=(H=|j-l| m=-j 

The free outgoing wave solutions to the vector Helmholtz equation are then krh^ (kr)Y? m (9, (f>). The vector spherical 
harmonics arc orthonormal in the usual way, 

[\m9d9 d4>Y^ mi {e,<t>y ■ Y£ m2 (9,<f>) = 6 jlj% S tlta S mim , (24) 
Jo Jo 

and under complex conjugation they transform as Y^ m (6,(j))* = (— iy+ e + m + 1 Yj_ m (9, (f). 

We can now use the basis of free spherical vector waves to set up the variable phase calculation in the same way 
as in the scalar case. In place of Eq. ©, we will need the integral over solid angle of the dot product of two vector 
spherical harmonics multiplied by a third ordinary spherical harmonic (since the potential is still expanded in terms 
of ordinary spherical harmonics), which is given in terms of the 6j-symbol and Clcbsch-Gordan coefficients as 

sinede [ * dcf>Y£ mi (9,<f>) ■Y£ ma (e,<i>)Yr(e,4>) = 

o Jo 



V Mi,_, V n /(2ji + l)(2j a + l)(2*i + l)(2& + l) \l x £ 2 £ 



In place of Eq. ([7j , we then have the coupling between channels 



, U i"+i'+e +m " +m ' + i (gj + l) (2j» + l) (2£+l) (2^ + 1) f ^ f rt 
"jll'ji'l" - \ l ) y 47r(2f + 1) \j" i ^ j 

With this modification, the calculation of the 5-matrix for the vector Helmholtz equation proceeds analogously to 
the scalar case. 



III. GENERALIZATION TO MAXWELL'S EQUATIONS 



To generalize to the case of electromagnetic scattering, we consider a linear, spatially-dependent dielectric with no 
free charge. The permittivity e(r) goes to one at large distances. We will treat each frequency to — cyW separately, 
so our formalism can easily incorporate frequency dependence in e(r), though as in the scalar case we do not indicate 
this possibility explicitly. The permittivity can also include an imaginary part, representing dissipation. We are 
interested in solutions to the Maxwell wave equation 

V x V x E k (r) = k\(r)E k {r) , (27) 

for k 7^ 0. Such solutions automatically obey Gauss's law V • Dk(r) = 0, where D k (r) = e(r)Ek(r). However, 
the solutions to this equation do not span the full space of vector functions, because in addition to these transverse 
solutions there also exist longitudinal solutions, which can be written as the gradient of a scalar function and therefore 
solve Eq. (|27[) with k — 0. This situation is problematic for the variable phase approach (in which we consider each 
k separately), because it implies that the matrix coefficient of the second derivative operator for fixed nonzero k will 
not be invertible, leading to an implicit differential-algebraic equation. We thus consider a modified equation that 
allows us to find the S-matrix for the transverse modes while avoiding this problem. 



A. Transverse and Longitudinal Modes 



To motivate our approach, we review a common method for solving the Maxwell wave equation in free space (or 
within a dielectric with constant permittivity), which is to replace the curl-curl operator V x Vx by minus the 
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Hclmholtz operator — V 2 . These operators commute, so they share the same eigenstates, and when acting on the 
transverse states, they share the same eigenvalues. (Recall that —V 2 E k (r) = VxVx E k (r) — V (V • E k (r)), where 
for transverse modes in empty space V • E k (r) = by Gauss's Law.) However, when acting on the longitudinal 
modes, the eigenvalue of — V 2 is the usual value of k 2 associated with a mode with wave number k, rather than zero. 
Once all the solutions to the Hclmholtz equation have been identified, it then is usually straightforward to discard 
the longitudinal modes and keep only the transverse modes. 

We now generalize this procedure for the case of a position-dependent dielectric. We first rewrite Eq. (|27l) in 
operator form as 

" VxVx...W(r) = fc 2 .E fe (r), (28) 



where . . . represents the argument of the operator. We then define the generalized Hclmholtz operator as 

1 



V X V X ... - V[V ■ (e(r) ■...)] E k (r) = k z E k {r) , (29) 
e(r) J 

which gives the same situation as in the free case: The operators in Eqs. (j28|) and (|29|) commute and share the 
same eigenstates. For the transverse modes, they share the same eigenvalues as well, but for the longitudinal modes, 
the eigenvalue of Eq. (|28|) is k 2 = 0, while the eigenvalue of Eq. (|29p is the usual nonzero value of k 2 associated 
with a mode of wave number k. We note that this approach would continue to work in the presence of a nontrivial 
permeability [i(r), with the only change being that V x Vx is replaced by V x -ttVx . 
We will thus solve for the S'-matrix associated with the wave equation 

V x V x E k (r) - e(r)V[V • (e(r) • E k (r))] = k 2 e(r)E k (r) . (30) 

Again, we decompose both the solution and the source in the appropriate spherical harmonic basis, 

oo j+l 



oo J + i 3 oo e 

Ek(r) = J2 E E -iWMll(M) <r) = J2 E Qw(r)Y/(M), 
j=o e=\j-i\ m=—j e'=om'=-£' 



(31) 

where Q m (r) goes to V^8^5 m Q at large r. As above, we denote the matrix outgoing wave solution, written in 
the vector spherical harmonic basis, by F k (r). We then substitute this expression into Eq. (|30|) and carry out the 
vector spherical harmonic algebra symbolically in Mathematica, using the identities in Appendix [A] to implement the 
differential operators and Eq. (|25p to carry out the convolution involved in multiplying by e(r). 
The result is an equation of the form 

-d 2 {k,r)F£(r)+d 1 {k,r)F' k (r)+d {k,r)F k {r)=Q, (32) 

where the matrices d$(k,r), d\{k,r), and d2(k,r) can depend on the dielectric profile and its derivatives, and prime 
denotes a derivative with respect to r. The replacement of Eq. (|28p by Eq. (f2Uf ensures that (^(fc, r) is an invertible 

matrix, so we let Di(k,r) = (d,2(k, r)\ di(k,r) and D (k 1 r) = (d2(k 1 r)^j do(k,r) to obtain 

-F£(r) + D 1 (k,r)F k -(r) + Do(k,r)F k (r) = 0. (33) 



Furthermore, since the generalized Hclmholtz operator in Eq. ((29]) approaches the ordinary Hclmholtz operator as 
e — > 1, for large r this equation approaches the ordinary Helmholtz equation, with Di(k,r) = and D (k,r) = 
— k ■ Again using Mathematica to carry out the symbolic algebra, we parameterize the outgoing solution by 
F k (r) = G k (r)W(kr) and, taking advantage of the simplifications arising from Eq. (fTTj) . obtain an ordinary matrix 
differential equation for G k (r), 



(r) + (D 1 (k,r)G k (r) - 2G' k {r)) (^\ogW(kr)\ +£>i(fc, r)G' k (r)+ (A>(fc, r) + fc 2 ) G k (r)-G k (r)^ =0, (34) 



with the boundary conditions G k (oo) = 1 and G" fc (oo) = 0. 

The solutions to Eq. (|30p include both the transverse solutions to the Maxwell equation that we are looking for and 
the longitudinal modes that wc wish to discard. Because the S'-matrix is defined in terms of incoming and outgoing 
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asymptotic waves, it is straightforward to project out the transverse modes. In the free case, the transverse solutions 
are given by [l9| 

M jm , k {r,6,4>) = Zj {kr)Y^ j (e,4>) (35) 



N jm , k (r, 0, 0) = -y^-^-iCfcr)^" 1 ^ <t>) + ]j ^~[Zj+i(kr)Y^ j+1 (0, 4>) (36) 

for j = 1,2,3..., where Zg(kr) is the appropriate spherical Bessel or Hankel function of order £. Since we have free 
electromagnetic waves far away from the dielectric, by simply projecting the S-matrix onto the subspace spanned by 
these transverse solutions at large distances, we obtain the full electromagnetic S-matrix. 



B. Inward/ Outward Integration in the Maxwell Case 



The presence of first-derivative terms in Eq. (|3"3")l necessitates some modifications of the Wronskian analysis that we 
used in Sec. Ill Bl to obtain the S'-matrix by combining the outgoing and regular solutions at an intermediate fitting 
point. We consider the transpose of the regular solution, obeying 

-*fe(r)'- ($ fc (r) t £> 1 (fc,r))' + $fc(r) t D (fc,r)=0. (37) 

which we again parameterize by $fc(r) t = W^r) -1 H k (r). We obtain the differential equation 

-Hk(r)+ (J^logWikrf) (# fc (r)£i (A, r) + 2H' k {r)) +2 (JLl og W(kr)j H k {r) 

-H' k (r)Di(k,r) - H k (r)D[{k,r) + H k (r) (po(k,r) + fc 2 ) -^H k ( r ) = 0, (38) 

with the boundary conditions H k (0) — and H' k (0) = 1. Now the quantity that is independent of r is not the 
Wronskian but instead 



W[$ k (r)\F k (r)}=W[$ k (r) t ,F k {r)}-$ k (r) t D 1 (k,r)F k {r). (39) 



Because the additional term in Eq. (|39[) vanishes at r = 0, the expression for the electromagnetic S'-matrix in terms 
of W is the same as in Eq. (fl9]l . with yV k \ r =r replaced by yV k \ r =r - 



IV. NUMERICAL RESULTS 



We have constructed "proof of concept" implementations of these calculations using Mathcmatica, which are avail- 
able from http://community.middlebury.edu/~ngraham. This high-level code provides a convenient illustration of 
our approach for small- to moderate-scale problems; more extensive calculations arc likely to require lower-level code 
making use of parallel linear algebra packages. In this section we describe sample calculations that use this code to 
verify and illustrate our approach. 



A. Consistency Checks 

Because some of the calculations we have described are the first of their kind, not all of our results can be compared 
with previous work. Nonetheless, we can verify a variety of complementary aspects of our calculations against known 
results or consistency conditions. In particular, we can check the following: 

• For potential scattering with real V(r) and electromagnetic scattering with real e(r), the S-matrix should be 
unitary, S|,S\. = 1, for real k. 

• For electromagnetic scattering, the S-matrix we obtain from solving Eq. (|30[) should commute with projection 
onto the asymptotic free transverse modes in Eq. (|36)) . 
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• For scalar, vector, and electromagnetic scattering, the result of the inward/outward calculation should be inde- 
pendent of the fitting point r . 

• For a spherical finite square well in the scalar case and a dielectric sphere in the electromagnetic case, the 
S'-matrix is diagonal and can be found analytically. For the scalar spherical square well, we have 0, p. 309] 

qhf ] (ka)f e (qa) - khf ] ' \ka)jt(qa) 

k ' £ ~ (Tj ny ' \ u > 

qh\ '(ka)j' e (qa) - kh\ (ka)jg(qa) 



where the potential is 



and q = \/k 2 + Vq . For the dielectric sphere, we have 0, p. 49] 

n 5 hf\ka)3z(nha) - hf^ (ka)j e (nka) 

ok,e,s = 7--7i) ttyt; , (.42; 

n°h\ (ka)j' e (nka) — h\ (ka)je.\nka) 

where we have defined the Riccati-Hankel functions ji(z) = zji(z), h^\z) = zh9 (z), and h^p(z) = zhg(z), 
5 = ±1 for the two transverse polarization channels, and the permittivity is 

n 2 r < a 



By using smooth functions that closely approximate the step functions in each case, we can verify that we obtain 
these results using our variable phase calculation. 

B. Sample Calculations 

To illustrate the numerical advantages of the variable phase method, we first consider a spherically symmetric 
example in electromagnetism, with 

/ n n~ c r f - , 1 — tanh [s(r — w)]\ .... 
ee m (r) = V^S m S m0 [l + h ^ ^ . (44) 



This profile gives a smooth approximation to a dielectric ball parameterized by height h, radius w, and edge steepness 
s. Because the profile is symmetric, the S'-matrix is diagonal and degenerate in the azimuthal quantum number 
to. Choosing our numerical matching point at ro = ^, we integrate outward starting from a small radius r smaU <C 

min {x,w) to obtain Hk{r) for r smaU < r < ro, and integrate inward starting from a large radius r big 3> max 
to obtain Gk(r) for r big > r > ro- Sample results are shown in Fig. [T] We see that these functions vary smoothly in 
response to the dielectric source, with trivial behavior outside the dielectric and no oscillations. In particular, Gfc(r) 
only becomes nontrivial when we reach values of r for which the source is no longer negligible; by choosing a moderate 
value of the steepness parameter s, we have softened the edge of the dielectric ball in order to highlight this transition. 

For comparison, we can reconstruct the normalized physical wavefunction il'T™ ( r ) f rom these results by writing 



G- 



k {r)W{-kr)P (W- k P - G k (r)W(kr)P (w k 

\ r—m \ 



C"( r ) = —?k= i v - r=ra - - v r=ro ' ( 45 ) 



P for r > ro 



W[kr)- x H k {r)C k for r < r 



where P is the projection matrix onto the transverse modes, the modified Wronskian Wk is evaluated at r = ro 
using Eqs. (|39p and (|16[) , and Ck is a constant matrix that matches the normalization of the two solutions, which is 
obtained by setting the two expressions in Eq. (|45[) equal at r = ro. This result, shown in Fig. [2j displays the typical 
oscillations associated with wave number k. By "factoring out" the free contribution W(kr), our method allows us 
to avoid these oscillations in numerical calculations. 
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FIG. 1: Eigenvalues of the matrices Gfc(r) and Hk(r) for k = 1, truncated at j max = 2, using the dielectric function in Eq. (|44|) 
with /i = 4, to = 1, and s = 8. For each eigenvalue, solid lines show the real part and dashed lines show the imaginary part. 
Taking ro = 5, we only calculate Hk{r) for r < ro and Gk{r) for r > ro- 




To illustrate the S-matrix as a function of fc, we consider a dielectric with a Drude model dependence on wave 
numb or, 

e tm (r) = V^S eo 5 m0 + - i==fi- - , ,„ Pim{r) , (46) 
— V-k 2 - (X p k) 2 

where p£ m (r) specifies the radial profile function for each spherical component of the dielectric profile. Here X p is the 
plasma wavelength, a v is the conductivity, and the frequency is uj = cVk 2 . We consider a deformed sphere using a 
profile given by 

r— 1 — tanh [sir — w) 1 , , 1 — tanh[s(r — w)] 
Pao {r) = V^tt p ^ and p 10 (r) = ^ ^ , (47) 

with all other pi m (r) equal to zero. The j = 1 eigenphase shifts for this case, given by one- half of the argument 
of the eigenvalues of the S'-matrix, are shown in Fig. [3] as functions of k. By comparing to the case where eoo(r) is 
kept the same but £io(0 is se t to zero, we see that a nontrivial eio(r) mixes the polarization channels and splits the 
degeneracy between \m\ = 1 and m = 0. As expected, these effects vanish at small k, where modes have wavelengths 
much larger than the length scale associated with the asymmetry, and also at large k, where modes have wavelengths 
much smaller than the plasma wavelength. 



V. DISCUSSION AND FUTURE DEVELOPMENTS 



We have developed a variable phase method to calculate the scattering S'-matrix for potentials in quantum mechanics 
and dielectrics in electromagnetism that are localized but do not have any particular symmetries. The result takes 
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FIG. 3: Eigenphase shifts, given by one-half the argument of the eigenvalues of the S-matrix, truncated at j max = 1. The left 
panel shows the case of the dielectric function given by Eqs. (|46(l and (|47(l . with A p = tt, a v = 1, w = 1, and s = 8, while the 
right panel shows the result for the same eoo(r), but with eio(^) = 0. 



the form of a matrix initial value ODE given in terms of a spherical harmonic decomposition of the scattering source. 
By using the Wronskian, we can combine inward and outward integration in r to obtain a well-behaved numerical 
computation, which remains tractable even for imaginary wave number k. Finally, we have extended this approach to 
the electromagnetic case by considering a modification of the Maxwell wave equation that avoids problems associated 
with disentangling the transverse and longitudinal waves. 

Our high-level Mathematica code provides a transparent and flexible high-level implementation of the methods 
described here, but it is only suitable for small- to moderate-scale calculations. Larger-scale calculations involving 
large numbers of partial waves will require the use of optimized low-level parallel linear algebra routines. Since the 
ultimate problem to be solved is quite generic, such calculations can take advantage of standard numerical packages 
for matrix ODEs. 
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Appendix A: Differential Operators 



Here we collect the differential operator relations needed to express Eq. (|30|) in the vector spherical harmonic 
basis, taken from Ref. [ill]. In these equations f(r) is an arbitrary radial function, YJ n (9, <f) is an ordinary spherical 



harmonic, and Y- m (6, cj>) is a vector spherical harmonic. 



V (/(r)17»(M)) 
V-(/(r)Y^ +1 (0,, 



1=3 ( 



Vx(/(r)Y^' +1 ( 



x (/(r)Y/= J '(M) 



3 + 1 



2 j + 1 \ dr r 
[7+1 ( d j + 2 



2j + 1 V* 



3-1 
2 j + 1 V dr r 

Id ,j + 2 



2j + 1 \dr 



3 



/(r)Y™(M) 

/(r)Y/"(M) 
f(r)Yff(9,<l>) 



3 + 1 f d j 



2j + 1 \ dr r 



3 



2j + 1 V dr r 



firm 



jm 1 









i ^ r J 



jm \ 
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V x 



{mYt?-\B,ti) = (| - ^) /(r)l£'(M) (Al) 



Appendix B: Free Green's Functions and Plane Wave Expansions 

Throughout this paper we have considered scattering in a spherical partial wave basis. For both Casimir calculations 
and traditional scattering problems, it is helpful to be able to convert these results to a plane wave basis. The key 
tools in this conversion are the expansion of a plane wave and the expansion of the free Green's function in terms of 
free spherical waves. Again drawing on Ref. JT8J , we collect those expansions here. For scalar scattering we have the 
well-known results 

e ik- r = ^^Mk^Yriek^kTYtie^) , (bi) 

lm 

where dk and <f>k are the angles of k in spherical coordinates, and 

Qo{ry,k)=ikY,Mto<)h?\kr > )Yr{6',<l/rxr{e,<l>), (B2) 

lm 

where r< (r>) is the smaller (larger) of r = \r\ and r' = \r'\. For vector waves, the expansion of a plane wave with 
polarization £ becomes 

€e*-'- = 47rX;i < (€.y J f n (fl fcl ^)*)^(Ar)^ ro (e,0), (B3) 

while the expansion of the free dyadic Green's function is 

G(r i , r 2 , k) = ik 3t (kr< (fcr> )Yf m (6 1 , & )* ® Yf m (0 a , fa ) . (B4) 

tjm 

We can also express these results in terms of transverse and longitudinal vector spherical harmonics. For the decom- 
position of a vector plane wave, we define 



for j = 1, 2,3 . . ., and 



where j = 0, 1, 2, 3 . . .. (Note that for j = 0, the unphysical term with £ = —1 is multiplied by zero.) Similarly, we 
consider the free transverse modes in Eqs. (|36j) along with the free longitudinal mode, given by 



L jm , k (r, M) = y ^jZi-iMY^- 1 ^, 4>) + y |^z J+1 (fcr)y/- J+1 ^, 0) (B7) 

for j = 0,l,2.... 

For the decomposition of a plane wave, we then have 

£e ik - = 4tt J2 ^° (€ ' X-^, fe (r, 0, 0) , (B8) 

Xj'm 

where cr = 0, 1, — 1 for x = M, A/", i respectively, and for the free dyadic Green's function we have 

G(n,r 2l k) = l kJ2 X%l k {r<T ® X$,*(r>) . (B9) 
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again for x = M, AT, L. Here the regular solution is given by taking zg(kr) = je(kr) in Eqs. (|36|) and (|B7|) . while the 
outgoing solution has Z[{kr) = h^{kr). 



[1] R. G. Newton, Scattering Theory of Waves and Particles (McGraw-Hill, New York, 1966). 

[2] K. Chadan and P. Sabatier, Inverse Problems in Quantum Scattering Theory (Springer- Verlag, Berlin, 1989). 

[3] E. I. Kats, Sov. Phys. JETP 46, 109 (1977). 

[4] M. T. Jaekel and S. Reynaud, J. Physique I 1, 1395 (1991). 

[5] T. Emig, N. Graham, R. L. Jaffe, and M. Kardar, Phys. Rev. Lett. 99, 170403 (2007). 

[6] T. Emig, N. Graham, R. L. Jaffe, and M. Kardar, Phys. Rev. D77, 025005 (2008). 

[7] S. J. Rahi, T. Emig, N. Graham, R. L. Jaffe, and M. Kardar, Phys. Rev. D80, 085021 (2009). 

[8] O. Kenneth and I. Klich, Phys. Rev. Lett. 97, 160401 (2006). 

[9] N. Graham, M. Quandt, and H. Weigel, Spectral Methods in Quantum Field Theory (Springer- Verlag, Berlin, 2009). 
[10] P. C. Waterman, Proceedings of the IEEE 53, 805 (1965). 
[11] P. C. Waterman, Phys. Rev. D 3, 825 (1971). 

[12] M. I. Mishchenkoa, G. Videenb, V. A. Babenkoc, N. G. Khlebtsovd, and T. Wriedte, J. of Quantitative Spectroscopy & 

Radiative Transfer 88, 357 (2004). 
[13] F. Calogero, Variable Phase Approach to Potential Scattering (Academic Press, New York, 1967). 
[14] M. T. H. Reid, A. W. Rodriguez, J. White, and S. G. Johnson, Phys. Rev. Lett. 103, 040401 (2009). 
[15] B. R. Johnson, J. Opt. Soc. Am. A 16, 845 (1999). 

[16] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C: The Art of Scientific 

Computing, 2nd edition (Cambridge University Press, Cambridge, 1992). 
[17] N. Graham, R. L. Jaffe, V. Khemani, M. Quandt, M. Scandurra, and H. Weigel, Nucl. Phys. B645, 49 (2002). 
[18] D. Varshalovich, A. Moskalev, and V. Khersonsky, Quantum Theory of Angular Momentum: Irreducible Tensors, Spherical 

Harmonics, Vector Coupling Coefficients, 3-j Symbols (World Scientific, Hackensack, NJ, 1988). 
[19] L. C. Biedenharn and J. D. Louck, Angular Momentum in Quantum Physics: Theory and Application (Encyclopedia of 

Mathematics and its Applications) (Cambridge University Press, Cambridge, 2009). 



